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, Abstract. Numerical schemes for Einstein's vacuum equation are developed. 

£N| ' Einstein's equation in harmonic gauge is second order symmetric hyperbolic. It is 

discretized in four-dimensional spacetime by Finite Differences, Finite Elements, and 
O ! Interior Penalty Discontinuous Galerkin methods, the latter related to Regge calculus. 

The schemes are split into space and time and new time-stepping schemes for wave 



£h . equations are derived. The methods are evaluated for linear and non-linear test 

problems of the Apples-with- Apples collection. 
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^ ! 1. Introduction 

■ Numerical methods for the solution of Einstein's equation in general relativity are 
mainly based on Finite Differences (FD) and Pseudo-Spectral-Collocation (Bonazzola 
et al. 2004, Boyle et al. 2007) schemes space so far. The Finite Element method (FEM), 
^ ■ or more generally Galerkin schemes have been used for reduced or auxiliary problems 
in numerical relativity (Arnold et al. 1998, Metzger 2004, Sopuerta et al. 2006, Aksoylu 
et al. 2008, Field et al. 2009). However, Galerkin methods are heavily used for the 
solution of wave problems in areas like acoustic and electro-magnetic scattering and 
elastic waves (Cohen 2002). This is mainly due to their way to deal with heterogeneous 
media and arbitrarily shaped geometric objects, represented by unstructured grids. 
Furthermore, the convergence theory of Galerkin methods is based on lower regularity 
(differentiability) requirements than Finite Differences and spectral methods. 

General relativity is governed by Einstein's equation, which can be written as a 
system of second order partial differential equations in spacetime. In order to define a 
well-posed initial-value (Cauchy) problem, additional gauge conditions are needed. For 
the numerical solution of the system, spacetime is usually split into space and time 3 + 1 
and finally a time-stepping scheme is derived. Using a lapse- and a shift-function, a 
sequence of space-like manifolds is constructed, which fixes the gauge freedom. There 
are many improvements of the original ADM (Arnowitt et al. 1962, York 1979) splitting 
like BSSN (Shibata & Nakamura 1995, Baumgarte & Shapiro 1999). The equations are 
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usually discretized in space by FD or spectral schemes and independently in time by an 
explicit integrator for ordinary differential equations. 

The harmonic approach and its generalizations first incorporate the harmonic 
gauge condition into Einstein's equation in spacetime to derive a hyperbolic system 
(Fock 1959, Bruhat 1962, Reula 1998, Friedrich k Rendall 2000, Pretorius 2005). 
Afterwards, the system is again split into space and time and discretized. Generalized 
harmonic methods modify the gauge condition, but usually preserve the hyperbolicity. 

In this paper, we follow a slightly different approach. Starting with the hyperbolic 
system of Einstein's equation in harmonic gauge, we discretize first in spacetime. 
Introducing a global time-step, the system is split afterwards in space and time. 
However, adaptive grid refinement in space and local time-stepping schemes can also be 
derived in a consistent way. This is similar to Regge calculus (Regge 1961, Sorkin 1975) 
in spacetime. 

The main contribution of the paper however is the development of a Finite Element 
and an Interior Penalty Discontinuous Galerkin (DG) method for Einstein's vacuum 
equation. Both methods are derived from a variational formulation, which is obtained 
from the Einstein-Hilbert action and harmonic gauge. In fact, Galerkin methods are 
always based on a variational version of the differential equations. 

Galerkin schemes have been considered for the discretization of wave equations in 
several ways so far: The wave equation duu = Au as an example problem is written in 
variational form as 

J n (d tt u)wd 3 x = -f n (Vu)-(Vw)d 3 xVw 

with trial functions w, integration over the spatial domain Q, and zero boundary 
conditions. This gives rise to FEM (Dupont 1973, Baker & Bramble 1979) and DG 
(Ainsworth et al. 2006, Grote et al. 2006, Field et al. 2009) in space schemes, used in 
conjunction with a standard time integrator like the leapfrog scheme. The first order in 
time formulation dtv = Au and dtu = v in variational version in time reads as 

— f T v(d t w)dt = J T (Au)wdtWw 

— f T u(d t w)dt = f T vw dt Ww 

on the interval T and without initial value terms. In order to obtain a time- 
stepping scheme, a time-discontinuous Galerkin method can be constructed (Jamet 
1978, Eriksson et al. 1985). Note that time continuous functions do not lead to a 
time-stepping scheme, but a single large equation system for all times. We can combine 
both Galerkin schemes to a spacetime FEM like 

J nxT v(d t w)dtd 3 x = f QxT (Vu) ■ (X7w)dtd 3 x Vw 

— f nxT u(d t w)dtd 3 x = f QxT vw dt d 3 x Ww , 

continuous (French & Peterson 1996, Anderson & Kimn 2007) and discontinuous 
(Hulbert & Hughes 1990, Monk & Richter 2005) in time. In this paper, however, we 
will consider second order in space and time formulations of type 

f nxT (d t u)(d t w)dtd 3 x = f nxT (Vu) ■ (Vw) dtd 3 x Vw , (1) 
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again without boundary and initial value terms. It can be re-written covariant and leads 
to time-stepping algorithms even for time-continuous Galerkin discretizations, which 
differ from first order formulations in general. 

The first result of the paper in section [2] is in fact the derivation of such a variational 
formulation of Einstein's equation from the Einstein-Hilbert action. In addition, a 
linearized formulation is discussed. 

If we restrict the solution and trial functions in ([1]) to some finite dimensional spaces, 
we obtain Galerkin discretizations in section [3j Although the spacetime formulation 
relates values at different points in space and time, it reduces to a time-stepping scheme 
for global time steps. The FEM scheme reduces further to the leapfrog time-stepping for 
piecewise linear functions in time, equidistant time steps, and without mixed space-time- 
derivatives. Note that leapfrog is related to the Stormer-Verlet scheme and a special 
case of the Newmark scheme. However, in the general spacetime case the FEM and the 
symmetric and non-symmetric DG spacetime schemes seem to be new. They form the 
next result of this paper, see sections 13.31 and 13.41 

While the leapfrog scheme is explicit for FD in space, see (Cohen 2002) and 
(Pretorius 2005, App. B), the FEM method in space requires the solution of a global 
equation system with mass matrix J n uw d 3 x each time step. The DG method in space is 
computationally more efficient than FEM in general, because the mass matrix is block- 
diagonal and the equation systems are easier to solve. However, by a special choice of 
numerical quadrature rules (mass-lumping) in FEM, see (Cohen 2002), and a choice of 
orthogonal ansatz functions in DG, see (Riviere 2008), the mass matrix is diagonal and 
the equation systems are trivial to solve. 

Now we put together the variational formulation of Einstein's equation and the 
spacetime Galerkin schemes and we obtain in section [331 as the main result, a FEM, a 
symmetric and a non-symmetric Interior Penalty DG method for Einstein's full vacuum 
equation. As an intermediate step we briefly discuss a simpler, linearized version of 
Einstein's equation. 

Memory requirements for nodal FD and piecewise linear FEM schemes for Einstein's 
equation are comparable, namely ten metric component values per grid node. The DG 
methods need this storage of 10 values for each element and each ansatz function, i.e. 
10 • 5 or 10 • 16 for linear or multi-linear functions, thus are more memory intensive. 
The fields are needed for two previous and the current time-slice in the leapfrog 
time-stepping. We put the discrete fields into the variational formulation, which now 
translates to non-linear equation systems. The matrix entries are computed by numerical 
quadrature rules. Additional storage may be required for the matrices and solution of 
the equation systems, which depends on the solver. 

Finally, some numerical experiments inspired by the Apples-with- Apples test suite 
(Alcubierre et al. 2004, Babiuc et al. 2008) are used to compare both schemes with a 
more traditional FD scheme in section |H The Galerkin schemes with piecewise linear 
functions on equidistant, cartesian grids show comparable CFL conditions, comparable 
second order accuracy, similar (sometimes opposite sign) dispersion second order in 
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grid spacing, and comparable second order accurate harmonic gauge conditions. The 
errors on unstructured grids additionally depend on the orientation of the elements with 
respect to the wave characteristics and element angle conditions. 

In order to solve realistic test cases in general relativity, techniques to handle 
apparent horizons are needed. Standard techniques include the puncture approach 
(Brandt & Bruegmann 1997, Campanelli et al. 2006, Baker et al. 2006), excision 
(Pretorius 2005), and singularity avoiding slicing conditions. Slicing would lead to 
a generalized harmonic gauge. Excision is compatible with harmonic gauge and the 
excised domain can be approximated by unstructured grids, which seems to be most 
promising. Furthermore, the Galerkin schemes have to be generalized to higher order, 
which is straightforward in space, but is more difficult in time for stability reasons. 



2. Einstein's Vacuum Equation 

2.1. Strong Formulation 

We start with the standard derivation of Einstein's equation via the Einstein-Hilbert 
action defined by 



s ■= / Ry/=<jdrx 

J M 

in the case of vacuum, in the notation of (Straumann 2004). We consider it as a function 
of the metric tensor g a p and its derivatives. The Ricci tensor R a p and the Ricci scalar 
R = g al3 R a {3 contain up to second order partial derivatives of g a p. We are looking for 
an extremum of S. The variation of S is 

5S= [ {R, v - L R )(5gnV=9d A x , (2) 

as long as the variation bg^ 11 vanishes at the boundary of the domain Ai. Otherwise we 
obtain an additional boundary term 

? f fsTOySgm - dp5 9lxu )n^d 3 x (3) 

z JdM 

which can be used later for boundary conditions using derivatives of g^ v . We rename 
the variation 

if ■= 8 g n» . 

The variational formulation reads as: We seek a solution g a p e V a such that 5S = for 
all v a p € V t with appropriate ansatz and trial spaces. Dirichlet boundary conditions on 
(parts of) dM. can be built into V a and V t : The solution takes the Dirichlet values and 
the trial functions vanishes there. Boundary conditions involving derivatives require an 
additional boundary term like ([3]). The variational formulation translates to the strong 
formulation as R^ — \g^, v R = or in vacuum 

Rfj,v = 
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with appropriate boundary conditions. However, in order to obtain a well posed initial- 
boundary-value or Cauchy problem, we need an additional gauge condition. We choose 
the standard harmonic gauge with 

r Q := <rr^ = o , (4) 

which is a condition on first order derivatives of g a p. This way, we can modify Einstein's 
equation as 

RjS ■= Rp, ~ ^aud^ - l -g a ,d u Y a = , (5) 
with principal part 

R$ PP ■= -\g aP d a d p g, u . (6) 

Now, we have a quasi-linear, second order, symmetric hyperbolic differential equation, 
which we will later discretize by finite differences. Note that this remains true if we 
switch to a generalized harmonic gauge. Equation (j4]) changes to T a = H a (x,g) with a 
gauge driver H. This driver may depend on coordinates and the metric, but must be 
independent of derivatives of g in order to preserve the principal part R^ pp . 



2.2. Variational Formulation 

Galerkin discretizations are based on a variational formulation. We start with the 
standard variational formulation ([2]). By Stokes' theorem, we can remove the second 
order derivatives. With harmonic gauge (jlj) we arrive at a variational version of ([HD 

a(g,v) := i f </°V=0 (d a ^)(c^)d 4 x , (7) 
M 



2 

which is symmetric in the first order derivatives of g^ and v ^ in the special case of a 
fixed background g^ u . Again there is an additional boundary term, if the variation v 
does not vanish on the boundary dAi 
1 

~ 2 

The remaining terms can be assembled in 

<l{9, v ) ■= \ I M 9 a ^9 pa V^9 ( {daQpn) (dpg^) - {d a g m ) (d a g Pu ) 

+ (d a g P ^(dugpa) + (du,gap)(d/3g*v) (9) 
-\ {d^g ap ){d u g^) ) v^d 4 x , 



g^y^(d^ u )v^d 3 x . (8) 

dM 



which is quadratic and symmetric in the first order derivatives of g^, compare also 
(Fock 1959, App. B). The variational formulation now reads as 

seek g 6 V a such that a(g, v) + q(g, v ) = Wv G V t , . 

and T a = . ^ ' 
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Note that metric g G V a in fllOp does not need to have well defined second derivatives 
as in (jSJ) and may be chosen in an appropriate Sobolev space. In the case of a non 
vanishing energy-momentum tensor T^ v additional terms of type 

Kd,v) := / (flwflW3 - l-g^ga^V^gT^ v^d 4 
Jm 1 

appear on the right-hand side of ( |T0i) . 

Different types of initial and boundary conditions can be imposed on dhA by 

standard procedures to define a Cauchy problem: Homogeneous Dirichlet values are 

directly incorporated into all functions in V a and V t . Inhomogeneous Dirichlet conditions 

are built into the solution g, either direct in the discrete numerical scheme, or via 

an additive splitting into a homogeneous auxiliary solution and a non-homogeneous 

function for the boundary conditions. Neumann boundary conditions and other 

conditions based on derivatives of the solution on parts of dAi lead to additional terms 

in a of type (jHJ), where dpg^ is replaced by the given derivatives. The functions in V a 

and Vt do not vanish there. "Natural" boundary conditions can be defined as vanishing 

term (|Sj), that is g al3 \f—g{dpg^ v )n£ > = 0. The conditions can be translated back into a 

strong formulation via (J3]). 

2.3. Linearized Equations 

In a weak field approximation of Einstein's equation, we neglect the first order derivatives 
in (El) and arrive some background metric g^ v '. In the variational version 

(TIP]) , we can neglect q(g, v) and solve for a(g, v) = instead, again for a fixed background 
metric g^ . 

a{g,v):=\f g<#>flj {d^^dpv^x (11) 
The linearized version of the harmonic gauge condition (j3J) reads 

g aP r u (d,g uP - l -d p g, v ) = . (12) 

Now, we simplify the equations even further and consider a weak field in flat space. 
The linearization is taken around Minkowski metric g = r\ := diag(— 1, 1, 1, 1) and we 
obtain the strong formulation 

- \u gia , = (13) 
with d a = 7] al3 dj3 and □ = d a d a . This translates to the variational version 

seek geV a such that a(g, v) := - / ^{d^^id^^x = \/v G V t .(14) 

2 Jm 

The harmonic gauge condition (fT2|) reduces to 

d^ u - ^V af3 d"g a p = , 
which can be further simplified by the substitution h^ v := g^ u — ^rj a,3 g al 3 to 

Wh^ = . (15) 
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The differential equation still is (ITS]) —^Dh^ = 0, now with a divergence free h. The 
gauge conditions are linear and can be incorporated into the spaces V a and V t . 

3. Numerical Schemes 

3.1. Finite Differences (FD) 

For illustration purposes, the first numerical spacetime scheme will be based on finite 
differences. We consider the discretization of a linear, scalar, second order wave 
equation — Du = with suitable initial and boundary conditions. On a one-dimensional, 
equidistant grid with grid spacing h, we choose the stencil h)— 2u(x)+u(x+h)) / h 2 , 

also abbreviated as [1 — 2 l]/h 2 , to approximate the second derivative. It is second 
order accurate for u smooth enough. The d'Alembert operator can be obtained by 
an application of the stencil along each coordinate axis on a cartesian grid. The two 
dimensional stencil at a grid point for example is 

Ui— \,j — + u i+l,j ~ ^ u i,j + u i,j+l 



h 2 h 2 
which gives the explicit time stepping scheme 



, 



using values at time slices i — 1 and i to calculate the values at time slice i + 1. This is 
the leapfrog scheme in time and can be written as 

u i+ x = 2ui - m-x + {h ) 2 A h Ui (16) 

with a FD approximation of the spatial derivatives A. Note that a CFL condition 
ho/hk < 1 for all k > must hold for stability reasons (Cohen 2002). The initial 
conditions can be prescribed at two times slices Xq = and Xq = h , the boundary 
values at Xfc = and x^ = 1. Modifications for other types of initial and boundary 
conditions do exist. 

3.2. Compact Finite Difference Stencils (FDM) 

In order to generalize the FD stencils to mixed first and second order derivatives, we 
consider an alternative construction. In the one dimensional case, first derivatives can 
be approximated by central stencils u'(x + h/2) w (u(x + h) — u(x))/h at grid points 
x + h/2. The second derivative can be calculated as a central stencil of first derivatives 
u"(x) ~ (u'(x + h/2) — u'(x — h/2))/h which reduces to the one- dimensional FD stencil. 
However, in two (and more) dimensions the construction differs, if we consider cell- 
centered first derivatives: We differentiate in one directions and average in the other 
direction(s): 

^ 1 ( u i-l,j - u ij , Ui-lj + x - IHj+1 
<A) M i+l/2. 1+1/2 ~ o i r 



hn hn 
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We obtain the second derivatives as stencils 

-1 2 -I 

! 

d d « 



(2h ) 



4 
2 



-2 
-1 



and <9 <9i 



4/i /ii 



-1 1 



1 -1 



The discretization of the d'Alembert operator again gives a time-stepping scheme for 
time slice i + 1. However, the scheme is no more explicit like (fT6l) . Let us write the 
difference stencil [1 2 l]/4 as the matrix M and the stencil [—1 2 — l]/(hi) 2 as matrix 
A. We obtain the scheme 

Mu i+1 = 2Mui - Mui_i - (h ) 2 Aui . (17) 

We can compute the values Ui + \ at time slice i + 1 by the solution of a linear equation 
system with matrix M using the values Ui and at time slices i and i — 1. The matrix 
is positive definite, symmetric, and of bounded condition number. Hence, the system 
is easy to solve numerically for large systems by standard iterative solvers. Again, the 
CFL condition limits the time step size ho- 



3.3. Finite Element and Petrov-Galerkin Methods (FEM) 

We start with a variational version of the d'Alembert operator (TjQ), a first step 
towards ffTTj) : 

I [ rf* '(a aM )(^)A = ovi) (is) 

Following standard procedures in FEM, we choose a set of global, continuous, piecewise 
polynomial ansatz and trial functions fa G V a and ipj G V% as a basis of finite dimensional 
spaces V a and Vj, and obtain a finite element method: Find the coefficients Ui of the 
solution u = J2iU l fa ^ V a , such that f|T8l) holds for all trial functions v G V t . This can 
be written in basis functions as 



and in matrix notation Au = with solution vector u and matrix A 



(19) 



Kl 



if* (9 a ^)(^Xi. 



.M 



This is a spacetime discretization. Introducing a global time step, we split functions 
fa[x) = 4>i(x )(f>j(xx, x 2 , x 3 ) and ^j, and the domain M. = T x Q into time and space. 
Further, mixed space-time derivatives r] a0 = rf b = do not occur with space index a, b. 
We obtain 

= \U T V Q \do^){do^)dt){^far 3 d"x) 
+ |(/ T 0^)(/ n ^(^)(^)rf 3 x) . 

We introduce the mass matrix M and matrix A by 

= U' n faip* d 3 x and 
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In order to solve a Cauchy problem with initial conditions, we deviate from standard 
FEM for self-adjoint problems in a single detail: In order to mimic the behavior of the 
spacetime FD schemes, we start with initial data at two time slices i — 1 and i and 
use the scheme to calculate the next time slice i + 1. We use piecewise linear functions 
<f>i(t) = max(l — \t — Xi\/ho, 0) and ip® in time for equidistant time-steps ho and obtain 
the system in time t^[ — 1 2 — 1}M +^[14 1]A, which is of leapfrog type 

(m - | A) u l+1 = (m + ^A^j Ul - (m - . (20) 

In the 1 + 1 spacetime case, piecewise linear functions on an equidistant space grid, we 
further obtain M = ^[1 4 1] and A = 2 - 1]. 

The method can be interpreted as a Petrov-Galerkin method with different ansatz 
V a and trial V t spaces: We use piecewise polynomial functions centered at a grid point 
i at time io and space location (ii, 12,^3) for a cartesian grid. The functions are chosen 
piecewise linear in time. Let the grid points be in the time domain (0, k) with initial 
conditions at iq = and io — 1. We compute the solution for all ansatz functions 
(f>i located at times (2,k). However, we choose the trial functions ipj located at times 
(1, k — 1). The trial functions lag behind one time slice, but are identical in space. This 
is exactly the idea to solve for the next time slice io + 1 and coincides with a leapfrog 
scheme for equidistant time steps. 

The solution of equation systems is the most expensive part of the time stepping 
procedure. The advantage of FD schemes for leapfrog (fl~6j) is that mass matrix M 
is the identity and no equation systems have to be solved. However, there is a 
common technique in FEM called "mass lumping" to obtain diagonal matrices M, too: 
Integration in terms of type f Q (pl^d^x and f T 0°0°c?t are approximated by numerical 
quadrature rules on each finite element. For piecewise (multi-) linear functions, 
quadrature rules prove to be sufficient, which are based on the function values at the 
element vertices only. This is the trapezoidal rule on an edge and its generalizations to 
rectangles, cubes, triangles and tetrahedra. The ansatz functions fulfill 4> S i(xj) = 5^ and 
the mixed products {(frffiMxk) = S^S jp. = for i ^ j vanish on all element vertices x^. 
Hence, off-diagonal entries TOy vanish and M is in fact a diagonal matrix. We arrive at 
the computational efficiency of an FD scheme f|T6l) . once the integration is done. 

Note that ( |T9i) defines spacetime FEM also for higher order methods in space by 
piecewise polynomial functions 4>l and ipj, by pseudo-spectral Galerkin schemes in space, 
on unstructured grids in spacetime, and for adaptive grid refinement in spacetime. The 
approach does not easily extend to higher order methods in time due to a lack of stability 
of the respective time-stepping schemes. 

3.4- Interior Penalty Discontinuous Galerkin Methods (DG) 

Again, we start with the variational problem ({TBI) . However, we choose piecewise 
polynomial ansatz and trial functions <pi and ipj, which are no longer continuous 
over element boundaries. This leads to additional terms. Consider a common face 
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Cij := dEiCidEj of two neighbor elements Ei and Ej and normal unit vector oriented 
from Ei to Ej. We denote the average by {u} := + («|^))/2 and the jump by 



it 



: (w|_bJ — on the face e^. Let the volume of the face be | e^- 1 . We split the 

integration over M. of 
grid. 

a(u, v) 



HI) into the integration over elements Ei and all faces of the 



(21) 



^2 ^« 



Cp 



-//" ; //;///'f J e ..MHfi 3 x = o 



The first jump term is obtained by Stokes' theorem, the second is added for reasons of 
symmetry of a, and the last term with penalty parameters c p and c e weakly imposes 
inter-element continuity. We have modified the penalty term, originally strictly positive 
for elliptic operators, by {f] Ti^n^g} due to the indefiniteness of the bi-linear form. 

We choose polynomial ansatz and trial functions on each element and combine 
them without continuity to global functions (pi and ipj. They define a basis of the finite 
dimensional spaces V a and V t . Find coefficients such that 

J2u l a^i,^ 3 ) = Vj. 

i 

The scheme is called the symmetric interior penalty discontinuous Galerkin scheme 
(SIPDG). Note that an opposite sign of the second jump term leads to the alternative 
non-symmetric NIPDG scheme, in our case with penalty c p = 0. Boundary conditions 
require modifications of the terms with outer boundary faces, see (Riviere 2008). 

If we use linear polynomials along each coordinate axis on an equidistant grid as 
before, we can calculate the difference stencils explicitly. In two dimensions n = 2 for 
example, we use the local nodal basis (1 — x )(l — £i), x (l — Xi), (1 — x )xi, x Xi and 
shift and scale it to each element. Again we solve for time slice i + 1 using slices i — 1 and 
i. However, now there are four degrees of freedom per element instead of one per node. 
With a penalty term c e 



1 and different constants c p in both directions, we obtain 



A X n = A* 
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and a 5-block scheme for the degrees of freedom in element at time i + 1 and position j 

AiflVt+ij = Ao-xUij-i + AofiUij + A ,xUi, j+ x - A-x,oUi-ij ■ (22) 

Note that for each element a linear equation system Ax a needs to be solved. It is of 
the size of number of ansatz functions, which is cheaper to solve than the single large 
equation system for the FEM. However, the amount of work can be further reduced: It 
is possible to choose the local ansatz functions orthogonal with respect to the bi-linear 
form such that Ax t o is in fact diagonal or even the identity and no systems need to be 
solved any more. This way, we obtain an explicit time-stepping scheme like (fT6j) . 

For a second order differential equation in time, we need two initial conditions, like 
u(0, Xx) and dou(0, Xx). This can be converted into data on two initial time slices i = 
and 101. However, for the DG schemes, we need an initial spacetime approximation in 
elements at times slices and 1. For a linear ansatz in time direction, initial data is 
needed at least at the beginning and end of both time slices, namely three initial values. 
These can be computed with a start-up calculation. 

3.5. Linearized Einstein's Equation 

In order to solve linearized Einstein's equation ffl3|) resp. ([HP , we can generalize the 
scalar schemes for Ou, apply these to each component g^, and set the background metric 
to Minkowski g = rj. The linear gauge condition ( Tl5jl needs to be fulfilled. Divergence- 
free initial data guarantees this for all times in the continuous case. However, numerical 
errors will lead to a violation of the gauge condition. DG methods easily allow for locally 
divergence-free ansatz functions on each element. In contrast, it is difficult to implement 
globally divergence-free symmetric tensor fields in FEM analogous to divergence-free 
vector fields for Maxwell's equation, see (Nedelec 1980, Nedelec 1986). 

In the case of a prescribed curved background metric, we have to solve the linear, 
variable coefficient problem R$ pp = 0. The FD stencils are no longer applicable and 
we switch to the compact FDM stencils. The FEM implementation is based on the 
variational formulation 

ij> 4 / g a ^\Rg {dMQM** = Vj • (23) 

1 i J M 

The DG method now reads as 

a(u,v) ■= lJ2iJ Ei 9 aP V-g (d a u)(d p v)d 4 x 
-I Ei<j /ey {9 a(3 V^9 <d p u} [v]d 3 x 

-lE^l^ig^V^gn^v^x [ } 

+ | J2i<j |^ Lj9 aP <4V=9}[u][v]d*x = 

where we have generalized the penalty term to {g n^rig^ —g}- The matrices M and A 
now depend on the background metric g af3 y/—g, which varies in spacetime. Procedures 
to construct a diagonal M like mass-lumping in FEM in section 13.31 and orthogonal 
ansatz functions in DG in in section 13.41 have to be performed on a per-element basis 
and are thus more expensive, as are procedures to construct divergence-free ansatz 
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spaces. Once the matrix entries have been computed, the linear equations system of 
type fl20|) and (122]) can be solved by standard solvers. 

3.6. Einstein's Vacuum Equation 

We generalize the compact FDM stencils to Einstein's vacuum equation ([5]): The 
variable metric g^ v and its second order derivatives R$ are chosen node centered (at 
grid points), but the first order derivatives are chosen cell centered. The inverse 
metric g^ u is used to calculate and T a and is also cell centered, defined as the inverse 
of the cell average of the metric g^ u . The products of averaged T enter the Ricci tensor, 
as well as the node centered derivatives of T. This way, we can use the standard formulas 
r « 7 := y^id^+d^-d^), R pv := 9.1;-^+^-^, (H, and © 
to set up non-linear, discrete Einstein's equation and derive the time-stepping scheme. 
Note that no code generated by a symbolic algebra program is needed. 

The FEM and DG Galerkin schemes can also be generalized to Einstein's equation. 
The form a ([7]) resp. (I2TI) and the quadratic term q ([9]) define the variational problem 
(flUk). The integration is done numerically. The integral j M is split into integrals over 
an element YIiIe ( an d a face in (l2Tj) ). The integrals over a single element Ei and 
face eij are approximated by a numerical quadrature rule. The integrands of a and q 
are evaluated at the quadrature points. 

The matrices M and A now depend on the current metric g af3 y/—g and the equation 
systems of type (I20p and (I2"2"j) are non-linear. The time-stepping schemes are implicit 
and require the solution of a non-linear equation system for each time-slice. The DG 
method leads to a set of easy to solve local equation systems for each element. The 
FDM and the FEM have global coupling of the degrees of freedom of a time slice. In 
both cases standard non-linear solvers can be used. Note that the explicit FD method 
both gives an initial guess for a locally fixed background metric g^ v and can be used as 
a preconditioner for the principle part in an iterative solver. 

The harmonic gauge condition (Tj0) now is a non-linear condition and cannot be 
incorporated into a linear ansatz space V a . Note that a change of variables leads to a 
formulation of Einstein's equation with a new metric := yf—gg^ and a linear gauge 
condition d^ u = 0, which could be built into V a . 

Note that Regge calculus also discretizes a variational principle in spacetime for 
simplicial grids (Sorkin 1975). It can be considered a DG spacetime scheme with piece- 
wise constant metric tensor g^ v . This way, fT2~TT) generalizes it to higher order and 
arbitrary element shapes. However, Regge calculus does not use coordinates and is 
based on purely geometric entities like edge lengths and defect angles. Furthermore, the 
variation is with respect to the degrees of freedom, which are the squared edge lengths 
in Regge calculus and values of the metric in ( 12TT) . 
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4. Applications 

4-1. Linear Plane Wave 

For illustration purposes, we perform some numerical experiments with the schemes of 
section^ The test cases are adapted from the Apples- with- Apples test suite (Alcubierre 
et al. 2004, Babiuc et al. 2008). We document and compare convergence and stability 
of the schemes in different settings. 

We start with a mono-chromatic traveling plane wave for linearized Einstein's 
equation (fT3l and ffT4"|) with harmonic gauge f|T5|) . We use periodic boundary conditions 
and a Courant factor 1/2. The one-dimensional (1+1) test case is defined on the 
spatial unit interval (0, 1(. The exact solution and initial data is goo — 9n — ~9oi — 
sin27r(xi— x ). We use an equidistant grid and run all schemes of sections l3~T1 to l3. 41 Note 
that the original Apples-with-Apples tests were constructed for non-linear numerical 
codes, such that very small wave amplitudes effectively ran the codes in the regime of 
the linearized equations. Standard non-linear solvers like Newton's method in this case 
reduce the problem to a linear one. Hence, we directly ran a linear code for the linear 
problem. This is why we can use arbitrary amplitudes of the solution rather than very 
small ones (Alcubierre et al. 2004). 

In figure [1] the time evolution of the spatial maximum error at the grid points for 
a resolution of h x = 1/200 is depicted for the FD, FDM, FEM, SIPDG and NIPDG. 
We use penalty parameters c p0 = 1 and c v \ = 2 for SIPDG. Note that continuous error 
norms like £2(0, 1) more natural for FEM show a similar behavior with exception of 
the very first time steps, where an additional interpolation error is added to the global 
error. The point-wise divergence is bounded, although we do not take any measures to 
control it. This does not seem to be necessary. The solution in figure [2] (left) shows the 
spatial errors at the final time xq = 1000. We see mainly dispersion and the phase error 
of the different schemes, no errors in the amplitude. This is why the error in fact even 
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linear wave in 1+1 dimensions 



linear wave in 2+1 dimensions 
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Figure 2. Solution of a linear plane wave in 1+1 for hi = 1/200 at final time xq = 1000 
(left). Evolution of the maximum error of a wave in 2+1 on a cartesian grid with 
hi = 1/100 (right). 
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Figure 3. Linear plane wave in 2+1 on some triangulations with hi = 1/100: evolution 
of the maximum error (left) and corresponding spatial grids (right). 



decreases after some time, see figure [T] (right) n = 25 and n = 50. We observe a second 
order convergence of the error, the phase error and the divergence in h\ for all schemes. 

The two-dimensional (2+1) test case is defined on the spatial unit square (0, 1( 2 
with periodic boundary conditions. The exact solution and initial data is goi = go2 = 
g 12 = sin27r(xi + x 2 - V2x 2 ), gu = g 2 2 = (v 7 ^ - l)g 01 , and #00 = y/2g 01 . We run 
all schemes on cartesian equidistant grids, see figure [2] (right), except for the NIPDG 
scheme for a lack of stability. The SIPDG penalty term is chosen as c e = 1/2, more 
precisely Cp , c = 1/Hq. The second order convergence is comparable to the 1 + 1 case. 

\ e ij\ e 

In order to test the dependence on the spatial grid, we run the FEM also on 
a number of triangular grids, both uniform (tri) and randomly distorted (tri*), see 
figure [31 Now we obtain a strong dependence of the error on the orientation of the 
elements. The longest element edges tangential to the direction of the wave leads to a 
larger approximation error than in normal direction or for quadratic elements. 
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Figure 4. Evolution of the maximum error of the linear robust stability test in 1+1 
for hi = 1/200 (left) and of a non-linear Gowdy wave in 1+1 for hi — 1/100 and 
hi = 1/200 (right). 

4-2. Robust Stability Test for Linear Waves 

Now we consider a stability test for the linear wave equation. The starting point is 
a random perturbation of the zero solution. We use periodic boundary conditions on 
(0, 1(, equal distributed [— e, e] random values for all initial data, with e := 2.5- 10~ 7 (/i3) 2 
according to (Alcubierre et al. 2004). In figure H] (left) we observe stability of all schemes 
with oscillatory solutions for NIPDG and compact stencil FDM. 



4-3. Nonlinear Polarized Waves in the Expanding Gowdy Universe 

The polarized Gowdy spacetime on the Torus T 3 is a model for a gravitational wave 
in an expanding universe (Gowdy 1971, New et al. 1998). We use periodic boundary 
conditions on the spatial unit interval (0, 1( in x% direction. The solution is constant 
along X\ and 22 direction. Since we use harmonic gauge, time axis Xq differs from 
(Alcubierre et al. 2004). We use a Courant factor 1/4. The solution g^ v is given by 

g = diag(-e (A+3a;o)/2 , e Xo+p , e x °~ p , e ^- Xo)/2 ) with 
p := Jo(27re z °)cos(27TX3) and 

A := -27re :co Jo(27re :co )Ji(27re :EO )cos 2 (27ra;3) - 27rJ (27r)J 1 (27r) 

+2(7re*°) 2 (J 2 ) (27re*°) + J?(27re*°)) - i(27r) 2 (J 2 (2vr) + J 2 (2vr)) 

We run schemes of section I3"To1 with 3rd order Gauss quadrature (two points in each 
coordinate direction) on an element. The SIPDG penalty terms are chosen as c p q = .5 
and c p i = 2. In figure H] (right) we see the error for spatial resolutions hi = 1/100 and 
hi = 1/200, which demonstrates second order convergence. The DG methods do not 
seem to be as stable as the others. However, many numerical schemes start to diverge 
at some time t due to the exponential growth of some of the solution components 
(Alcubierre et al. 2004). 



FEM, DG, and FD Evolution Schemes in Spacetime 



16 



Conclusion 

We have developed new spacetime Finite Element (FEM) and Interior Penalty 
Discontinuous Galerkin (SIPDG and NIPDG) schemes for second order symmetric 
hyperbolic wave equations. The Discontinuous Galerkin schemes are computationally 
more efficient, but require more memory than FEM and Finite Differences methods. 
A variational formulation of Einstein's equation in harmonic gauge was derived, based 
on up to first derivatives of solution and trial functions. This led to new Galerkin 
schemes for numerical relativity. The schemes were presented and tested for second 
order accurate Galerkin schemes with multi-linear functions and global time steps. The 
Gowdy wave test demonstrated the need for additional numerical stabilization. This 
might be obtained by spatial filtering, artificial viscosity, or streamline diffusion. 

Extensions to arbitrary spacetime grids or (adaptive) local grid refinement in 
spacetime are straightforward, but may lead to larger and more expensive to solve 
equation systems. Higher order polynomials or other more accurate (spectral) function 
spaces improve the spatial accuracy of the schemes. However higher order in time 
schemes are more difficult to construct. 
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